/**
  ETFAna project, Anyang Normal University && IMP-CAS
  \class ETFVDC
  \brief virtual class for DCs in general - just for abstraction
  \author SUN Yazhou, asia.rabbit@163.com
  \since 2021-11-01
  \date 2024-03-07 last modified
  \attention
  changelog
  <table>
  <tr>  <th>Date         <th>Author      <th>Description                   </tr>
  <tr>  <td>2021-11-01   <td>Asia Sun    <td>file created                  </tr>
  </table>

  \copyright Copyright (c) 2021-2024 Anyang Normal U. && IMP-CAS with LGPLv3 LICENSE
*/

#include <TClonesArray.h>
#include <TH1F.h>
#include <TGraph.h>
#include "ETFVDC.h"
#include "ETFDetectorManager.h"
#include "ETFHPTDCCh.h"
#include "ETFVTrack.h"
#include "ETFMath.h"
#include "ETFParaManager.h"
#include "ETFMsg.h"
#include "ETFDCHit.h"
#include "ETFYaml.h"

#define er ETFMsg::Error

ETFVDC::ETFVDC(const char *name) : ETFDetector(name), type{0}, fNL(-1),
    fSenseWireDistance(-999.), fMaxDL(-999.), fLayerSpace12(-999.),
    fLayerSpaceXY(-999.), fStagger12(-999.), fX0Min(-999.),
    fDTLB(-999.), fDTHB(-999.), fhrt(0){
  const ETFYaml &str = (*fpm->GetMisPara())["STR"];
  if(str && str[name]) fvstr = str[name].as<vector<double>>();
  SetDCId(atoi(name + strlen(name) - 1));
  fhrt = fpm->GetRTHist(ETFMsg::Form("hrt_%s", name));
} // end ctor

ETFVDC::~ETFVDC(){}

int ETFVDC::Multi(DCType typeDC, int x12){
  return GetHitsArray(typeDC, x12)->GetEntriesFast();
} // end member function Multi

// add the data channel to the detector
void ETFVDC::Assign(const ETFRawChannel *ch, unsigned uid){
  const int wireId = fNL*(uid>>10 & 0x1) + (uid>>11 & 0xFF);
  if(!fIsFired) fIsFired = true;
  AddHit(ch, wireId, DCType(uid>>8 & 0x3));
} // end member function Assign

// generate uid
unsigned ETFVDC::uid(unsigned typeDC, int x12, int wireId){
  if((x12 != 0 && x12 != 1) || wireId < 0 || wireId >= fNL)
    er("ETFVDC", "uid: invalid input.");
  return fDetId+(typeDC<<8)+(x12<<10)+(wireId<<11);
} // end member function uid

// add DC signal wires to fGDumb
void ETFVDC::ConstructGraph(){
  if(fGDumb->GetN()) return; // has called

  fGDumb->SetMarkerStyle(7); fGDumb->SetMarkerColor(15); // gray
  for(int i = 0; i < 2; i++) for(int j = 0; j < fNL; j++){
    double p[3] = {fSenseWireDistance*j+fStagger12*i+fX0Min, 0., fLayerSpace12*i};
    TransformToLabRef(p);
    fGDumb->SetPoint(i*fNL+j, p[2], p[0]);
  } // end for over i
} // end member function ConstructGraph

/// \retval false: if obtained drift time out of range
/// NOTE that if trk has non-FitStatus::kUNFITTED FitStatus, this method
/// would do nothing, as drift distance only changes with drift time,
/// whose changes would ensure fFitStatus be set to FitStatus::kUNFITTED
bool ETFVDC::SetDriftDistance(ETFVTrack *trk){
  // drift distance already assigned
  if(trk->GetFitStatus() != FitStatus::kUNFITTED) return true;
  if(-999. == fDTLB || -999. == fDTHB)
    er("ETFVDC", "SetDriftDistance: drift time bound not set yet");

  for(auto &p : trk->GetHits()) if(!SetDriftDistance(p.second)) return false;

  return true;
} // end member function SetDriftDistance

bool ETFVDC::SetDriftDistance(ETFDCHit *h){
  if(!ETFMath::Within(h->GetDriftTime(), fDTLB, fDTHB)) return false; // false hit
  h->SetDriftDistance(SpaceTimeRelation(h->GetDriftTime()));
  return true;
} // end member function SetDriftDistance

// r = r(t); drift distance - drift time relation
double ETFVDC::SpaceTimeRelation(double t){
  if(t < 0.) return 0.;
  if(fhrt){
    if(t > fDTHB) return fMaxDL;
    return fhrt->Interpolate(t);
  } // end if
  return ParaSTR(fMaxDL, fDTHB*0.6, t);
} // end member function SpaceTimeRelation

// initial space-time relation (STR) using parabola passing (0,0), (rm,tm)
double ETFVDC::ParaSTR(double rm, double tm, double t){
  if(rm < 0. || tm < 0.) er("ETFVDC", "ParaSTR: minus rm=%f, tm=%f", rm, tm);

  const double a = -rm/(tm*tm);
  const double b = 2.*rm/tm;
  const double r = (a*t+b)*t;
  return r;
} // end member function ParaSTR

// finish the processing of all the assigned data
void ETFVDC::PostAssign(){
  ETFDetector::PostAssign();
} // end member function PostAssign


/////////////////////// THE CONTROVERSIAL FILTER METHOD ////////////////////////
// only 1 wire of adjacent wires is supposed to be fired: this method filter the others
// it is now found that for the genuine hit the tl is earlier and tot is much larger
// compared with the secondary spurious hits
// for a continuous range of wires, only 1 will be selected //
#define CONSERVATIVE_FILTER
#ifdef CONSERVATIVE_FILTER
// -- Mar 7, 2024: we only select those with earliest time and largest tot at the
// same time, or we don't make any selection at all (i.e., we let'em all pass)
void ETFVDC::Filter(TClonesArray *c){
  if(fPos[2] < 0.) return; // do not apply this to DCs around TA - just does not work
  const int n = c->GetEntriesFast();
  if(n < 4) return; // no need for any selection
  double tot0 = -999., totm = -999.; // the tot estimator of the current and optimal wire
  double t0 = -999., tm = -999.; // the time estimator of the current and optimal wire
  int i0 = 0; // the start of a continusous sub-range
  ETFDCHit *h0 = nullptr, *h1 = nullptr; // of the last and current hit
  ETFDCHit *hmtot = nullptr, *hmt = nullptr; // optimal hit of tot and time
  for(int i = 0; i < n; i++){
    h1 = dynamic_cast<ETFDCHit *>(c->At(i));
    if(!h1) er("ETFVDC", "Filter: not a DC hit");
    const ETFHPTDCCh *t = h1->GetTChannel();
    if(!t) continue;
    tot0 = t->tot();
    if(-999. == tot0) tot0 = 5.; // maybe too small to be measurable
    t0 = t->time();
    if(-999. == t0) er("ETFVDC", "Filter: abnormal timing signal for %s", GetName());
    t0 = 1000000. - t0; // so as to be positive, and the earlier the larger
    if(0 == i){ // the 1st hit
      totm = tot0; tm = t0;
      hmt = hmtot = h0 = h1;
      continue;
    } // end if
    // ATTENTION: FIXME FIXME FIXME: THIS FUNCTION RELIES ON HITS STORED CONTINUOUSLY
    // IN TERMS OF WIER ID IN THE TCLONESARRAY
    if(abs(h0->GetWireId() - h1->GetWireId()) == 1){ // adjacent wires
      if(t0 > tm){ hmt = h1; tm = t0; }
      if(tot0 > totm){ hmtot = h1; totm = tot0; }
    } // end if
    else{ // continuity broken
      // see what we've got //
      if(hmtot == hmt && i - i0 > 1){
        for(int j = i0; j < i; j++) dynamic_cast<ETFDCHit *>(c->At(j))->SetFiredStatus(0);
        hmt->SetFiredStatus(1);
      } // end if
      // and then start a new continuous range
      i0 = i;
      totm = tot0; tm = t0;
      hmt = hmtot = h1;
    } // end else
    h0 = h1;
  } // end for over DC hits

  // see what we've got //
  if(hmtot == hmt && n - i0 > 1){
    for(int j = i0; j < n; j++) dynamic_cast<ETFDCHit *>(c->At(j))->SetFiredStatus(0);
    hmt->SetFiredStatus(1);
  } // end if
} // end member function Filter

#else

// ---- the old one - Jan13, 2024 ---- //
// only 1 hit will be left in a range of continuous hits in all situations
// the estimator: arrive early + 0.005*tot //
void ETFVDC::Filter(TClonesArray *c){
  if(fPos[2] < 0.) return; // do not apply this to DCs around TA - just does not work - Mar 7, 2024
  double q1 = -999., qm = -999.; // the estimator of the current and optimal wire
  ETFDCHit *h0 = nullptr, *h1 = nullptr, *hm = nullptr; // of the last, current and optimal wire
  for(TObject * b : *c){
    h1 = dynamic_cast<ETFDCHit *>(b);
    if(!h1) ETFMsg::Error("ETFVDC", "Filter: not a DC hit");
    // if(!h1->fired()) continue;
    const ETFHPTDCCh *t = h1->GetTChannel();
    if(!t) continue;
    double tt = t->tot();
    if(-999. == tt) tt = 5.; // maybe too small to be measurable
    q1 = tt*0.005 - t->time(); // the larger, the better
    if(-999. == qm){ qm = q1; h0 = hm = h1; continue; }
    // ATTENTION: FIXME FIXME FIXME: THIS FUNCTION RELIES ON HITS STORED CONTINUOUSLY
    // IN TERMS OF WIER ID IN THE TCLONESARRAY
    if(fabs(h0->GetWireId() - h1->GetWireId()) == 1){ // adjacent wires
      if(qm < q1){
        hm->SetFiredStatus(0); // nullify h0
        qm = q1; hm = h1;
      } // end if
      else h1->SetFiredStatus(0); // nullify h1
    } // end if
    else qm = -999.; // start a new continuous wire range
    h0 = h1;
  } // end for over t0arr strips
} // end member function Filter
#endif
